'''
Created on Jul 22, 2009

@author: mkiyer
'''

from veggie.align.alignment import Alignment
import sys

class BowtieParser(object):
    '''
    The bowtie aligner outputs each alignment on a separate line. Each line 
    is a collection of 8 fields separated by tabs; from left to right, 
    the fields are:

    1. Name of read that aligned
    2. Orientation of read in the alignment, - for reverse complement, 
        + otherwise
    3. Name of reference sequence where alignment occurs, or ordinal ID if 
        no name was provided
    4. 1-based offset into the forward reference strand where leftmost 
        character of the alignment occurs
    5. Read sequence (reverse-complemented if orientation is -)
    6. ASCII-encoded read qualities (reversed if orientation is -). 
       The encoded quality values are on the Phred scale and the 
       encoding is ASCII-offset by 33 (ASCII char !).
    7. Number of other instances where the same read aligns against the same
       reference characters as were aligned against in this alignment. This 
       is not the number of other places the read aligns with the same 
       number of mismatches. The number in this column is generally not a 
       good proxy for that number (e.g., the number in this column may be 
       '0' while the number of other alignments with the same number of 
       mismatches might be large). This column was previously described as 
       "Reserved".
    8. Comma-separated list of mismatch descriptors. If there are no 
       mismatches in the alignment, this field is empty. A single descriptor 
       has the format offset:reference-base>read-base. The offset is 
       expressed as a 0-based offset from the high-quality (5') end of 
       the read.
    '''
    def __init__(self):
        pass
    
    def parse_file(self, thisfhd):
        i = 0
        m = 0
        prev_name = ''
        mappings = []

        for thisline in thisfhd:
            # count lines for status updates
            i+=1
            if i == 100000:
                m += 1
                sys.stderr.write(" %d\n" % (m*100000))
                i=0
            # check to make sure line can be parsed
            if not thisline:
                raise ParseError("line is None")
            thisline = thisline.strip()
            if not thisline:
                raise ParseError("blank line")
            if thisline.startswith("#"):
                # comment line
                continue
            
            # split the line into fields and extract relevant information
            fields = thisline.split('\t')
            try:
                name = fields[0]
                strand = fields[1]
                # convert the chromosome mapping field to ucsc chromosome format
                chrom = fields[2].split('|')[-1]
                if chrom.startswith('NC'):
                    chrom = str(int(chrom.split('_')[-1]))                    
                    if chrom == '23':
                        chrom = 'X'
                    elif chrom == '24':
                        chrom = 'Y'
                    chrom = 'chr' + str(chrom)
                chrom_start = int(fields[3]) - 1
                chrom_end = chrom_start + len(fields[4])
                # TODO: handle mismatches
            except:
                print thisline
                sys.exit(1)

            # alignment
            if chrom:
                if len(mappings) > 0 and name != prev_name:
                    yield mappings
                    mappings = []
                prev_name = name
                mappings.append(Alignment(chrom, chrom_start, chrom_end, strand))
        if len(mappings) > 0:
            yield mappings
                
